function [a_sims_all,cstar_sims_all,apstar_sims_all, age_sims_all, mpc_sims_all, mpc_sims_all2, mpc_sims_all3, mpc_sims_all4] = ...
    bc_tremble_sim(w,param,y_sims_all,reset_shocks,eps_eta_mat,sigma_tremble, b,cstar,asims0)

% INPUT
%   param:          parameter vector
%                   [bta,r,abar,W_cc,kappa,delta,deltax,ybar,sigma_y,
%                   cstar_ya,cbar,prior_mean_coeff,sigma_c,psi,cov_fun]
%   y_sims_all:     simulated y                             [sim_periods,nsim]
%   reset_shocks:   reset shocks                            [sim_periods,nsim]
%   b:              borrowing constraint (<=0)              [1,1]

% OUTPUT
%   a_sims_all:             simulated values for a              [sim_periods,nsim]
%   chat_sims_all:          simulated values for chat           [sim_periods,nsim]
%   cstar_sims_all:         simulated values for cstar          [sim_periods,nsim]
%   aphat_sims_all:         simulated values for aphat          [sim_periods,nsim]
%   apstar_sims_all:        simulated values for apstar         [sim_periods,nsim]
%   eta_sims_all:           simultaed values for eta            [sim_periods,nsim]
%   sigma_etasq_sims_all:   simulated values for sigma_etasq    [sim_periods,nsim]

% unpack parameters
bta = param(1);
r = param(2);
abar = param(3);
W_cc = param(4);
kappa = param(5);
deltax = param(6);
sigma_c = param(7);
psi = param(8);
cov_fun = param(9);
gam = param(10);
sigma_l = param(11);
lbar = param(12);
delta_eta = param(13);
delta_inh= param(14);
factor_mpc1 = param(15); 
factor_mpc2 = param(16); 


[sim_periods, nsim] = size(y_sims_all);
%{
if lognormal == 1
    [xs,~,cstar,~,ctilde] = bc_opt_pol_iid_egm_logn(bta,r,b,ybar_logn,(sigma_y^2+sigma_yi^2)^0.5,gam,amin,amax,na, util_type);
else
    [xs,~,cstar, ~, ctilde]    = bc_opt_pol_iid_egm(bta,r,b,ybar     ,(sigma_y^2+sigma_yi^2)^0.5,gam,amin,amax,na, util_type);
end


% set prior
%[xs,apstar,cstar] = bc_opt_pol_iid_egm(bta,r,b,ybar,(sigma_y^2+sigma_yi^2)^0.5,gam,amin,amax,na, util_type);
% this is to keep cstar below the balanced consumption level;
[xs,xi] = unique(xs);
temp = griddedInterpolant(xs,cstar(xi),'linear','linear');
%cstar = @(x) temp(x);

stateGrid = linspace( xs(1), xs(end),10001)';
stateGrid_step = stateGrid(2) - stateGrid(1);
cstar_interp = temp(stateGrid);
cstar = @(x) fast_interp( x, cstar_interp, stateGrid_step, stateGrid,10001);
%}

a_sims_all          = NaN(sim_periods,nsim);
cstar_sims_all      = NaN(sim_periods,nsim);
apstar_sims_all     = NaN(sim_periods,nsim);
age_sims_all        = NaN(sim_periods,nsim);
mpc_sims_all        = NaN(sim_periods,nsim); 
mpc_sims_all2        = NaN(sim_periods,nsim); 
mpc_sims_all3        = NaN(sim_periods,nsim); 
mpc_sims_all4        = NaN(sim_periods,nsim); 

for j=1:nsim
    %  preallocate
    y_sims               = y_sims_all(:,j);
    a_sims               = NaN(sim_periods,1);
    cstar_sims           = NaN(sim_periods,1);
    apstar_sims          = NaN(sim_periods,1);
    age_sims             = NaN(sim_periods,1);
    mpc_sims             = NaN(sim_periods,1); 
    mpc_sims2            = NaN(sim_periods,1); 
    mpc_sims3            = NaN(sim_periods,1); 
    mpc_sims4            = NaN(sim_periods,1); 
    eps_tremble_sims     = eps_eta_mat(:,j); 

    %  initial states
    a_sims(1) = asims0(j);
    age_sims(1) = 1;
    
    cstar_sims(1)    = min(max(cstar(a_sims(1)*(1+r) + w*y_sims(1)) + sigma_tremble*eps_tremble_sims(1),1e-3), a_sims(1)*(1+r) + w*y_sims(1) - b);
    apstar_sims(1) = w*y_sims(1)+ (1+r)*a_sims(1) - cstar_sims(1);
    
    %  period 2+
    for i = 2:sim_periods
        if reset_shocks(i, j) < deltax
            a_sims(i)   =  (w*y_sims(i-1) + (1+r)*a_sims(i-1) - cstar_sims(i-1))*delta_inh;
            
            age_sims(i) = 1;
        else
            a_sims(i)    = w*y_sims(i-1) + (1+r)*a_sims(i-1) - cstar_sims(i-1);
            age_sims(i)=age_sims(i-1)+1;
            
        end
        
        cstar_sims(i)    = min(max(cstar(a_sims(i)*(1+r) + w*y_sims(i)) + sigma_tremble*eps_tremble_sims(i),1e-3) , a_sims(i)*(1+r) + w*y_sims(i) - b);
        apstar_sims(i) = w*y_sims(i) + (1+r)*a_sims(i)-cstar_sims(i);
        
        y_mpc=y_sims(i)*(1 + factor_mpc1);
        cstar_mpc    = min(max(cstar(a_sims(i)*(1+r) + w*y_mpc) + sigma_tremble*eps_tremble_sims(i),1e-3), a_sims(i)*(1+r) + w*y_mpc - b);
        mpc_sims(i)=(cstar_mpc-cstar_sims(i))/(w*factor_mpc1*y_sims(i));                                   
        
        y_mpc=y_sims(i)*(1 - factor_mpc1);
        cstar_mpc    = min(max(cstar(a_sims(i)*(1+r) + w*y_mpc) + sigma_tremble*eps_tremble_sims(i),1e-3), a_sims(i)*(1+r) + w*y_mpc - b);
        mpc_sims2(i)=(cstar_mpc-cstar_sims(i))/(w*factor_mpc1*y_sims(i));                                        

        y_mpc=y_sims(i)*(1 + factor_mpc2);
        cstar_mpc    = min(max(cstar(a_sims(i)*(1+r) + w*y_mpc) + sigma_tremble*eps_tremble_sims(i),1e-3), a_sims(i)*(1+r) + w*y_mpc - b);
        mpc_sims3(i)=(cstar_mpc-cstar_sims(i))/(w*factor_mpc2*y_sims(i));   
        
        y_mpc=y_sims(i)*(1 - factor_mpc2);
        cstar_mpc    = min(max(cstar(a_sims(i)*(1+r) + w*y_mpc) + sigma_tremble*eps_tremble_sims(i),1e-3), a_sims(i)*(1+r) + w*y_mpc - b);
        mpc_sims4(i)=(cstar_mpc-cstar_sims(i))/(w*factor_mpc2*y_sims(i));   
            
    end
    
    a_sims_all(:,j)         = a_sims;
    cstar_sims_all(:,j)     = cstar_sims;
    apstar_sims_all(:,j)    = apstar_sims;
    age_sims_all(:,j)= age_sims;
    mpc_sims_all(:,j)=mpc_sims;
    mpc_sims_all2(:,j)=mpc_sims2;
    mpc_sims_all3(:,j)=mpc_sims3;
    mpc_sims_all4(:,j)=mpc_sims4;

end

end